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Estimating  damage  in  structural  systems  is  a  challenging  problem  due  to  the  complexity  of  the  likelihood 
function  describing  the  observed  data.  From  a  Bayesian  perspective  a  complicated  likelihood  means  effi¬ 
cient  sampling  of  the  posterior  distribution  is  difficult  and  standard  Markov  Chain  Monte  Carlo  samplers 
may  no  longer  be  sufficient.  This  work  describes  a  population-based  Markov  Chain  Monte  Carlo  approach 
for  efficient  sampling  of  the  damage  parameter  posterior  distributions.  The  approach  is  shown  to  accu¬ 
rately  estimate  the  state  of  damage  in  a  cracked  plate  structure  using  simulated,  free-decay  response 
data.  The  use  of  this  approach  in  identifying  structural  damage  has  not  previously  been  explored. 
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1.  Introduction 

The  problem  of  identifying  damage  in  structural  systems  can  be 
effectively  cast  as  an  estimation  problem.  That  is  to  say:  given  a  se¬ 
quence  of  observed  data  collected  from  one  or  more  sensors,  esti¬ 
mate  the  damage  presence,  location,  and  extent.  This  information 
is  necessary  for  predicting  the  future  state  of  the  structure  in  ques¬ 
tion,  i.e.  for  prognostics. 

There  are  two  main  classes  of  estimation  approaches  available: 
the  method  of  maximum  likelihood  and  Bayesian  estimation.  The 
former  produces  a  single  “best”  estimate,  defined  as  that  which 
maximizes  the  probability  of  having  observed  a  given  set  of  data 
(e.g.  structures’  vibrational  response).  The  latter  views  the  param¬ 
eter  to  be  estimated  as  a  random  variable  and  the  goal  is  to  com¬ 
bine  the  likelihood  with  prior  information  to  produce  the 
estimated  probability  density  function  (PDF)  associated  with  each 
parameter. 

The  difficulty  with  either  approach  in  structural  estimation 
problems  is  that  the  likelihood  function  (the  core  of  both  estima¬ 
tion  approaches)  is  often  an  extremely  complicated  function  of 
the  structure’s  parameters  with  many  maxima.  This  severely 
restricts  the  use  of  standard  optimization  algorithms.  For  example, 
a  simple  gradient  ascent  algorithm  cannot  be  used  to  find  the  max¬ 
imum  likelihood  estimate  (MLE)  in  many  structural  system  identi¬ 
fication  problems.  This  has  been  recognized  by  a  number  of 
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researchers  and  solutions  proposed.  Horibe,  for  example,  used  a 
genetic  algorithm  (GA)  to  find  the  MLE  of  a  structure’s  parameters 
in  a  simple  vibration  problem  [1  ].  Likewise,  Stull  et  al.  considered  a 
GA  in  generating  MLEs  of  structural  parameters  in  a  static  shell 
buckling  problem  [2].  In  the  context  of  damage  detection  Panigrahi 
et  al.  used  a  GA  to  perform  parameter  estimation  by  minimizing 
differences  between  predicted  and  observed  modal  properties  [3] 
using  a  likelihood  based  on  modal  properties.  Additional  work  by 
Hwang  et  al.  also  employed  a  cost  function  based  on  modal  prop¬ 
erties  and  used  a  GA  to  estimate  stiffness  in  a  composite  structure 
[4].  The  authors  point  out  that  for  certain  problems  maximizing  the 
likelihood  does  not  necessarily  require  a  GA.  For  example,  if  en¬ 
ough  is  known  about  the  parameter  distributions  a  priori  and  one 
has  a  good  initial  guess,  the  perturbation  approach  of  Fonseca 
et  al.  [5]  or  Xu  et  al.  [6]  can  be  used  to  obtain  MLEs  of  a  structure’s 
parameters.  In  the  related  field  of  structural  reliability,  researchers 
have  also  long  recognized  the  challenges  in  structural  optimization 
problems.  The  early  work  of  Kiureghian  and  Dakessian  [7],  for 
example,  developed  a  method  for  avoiding  local  optima  in  estimat¬ 
ing  failure  probability  in  structures.  More  recently  Guo  et  al.  [8] 
developed  an  algorithm  for  efficiently  obtaining  globally  optimal 
solutions  for  structural  design  problems  possessing  multiple  local 
optima. 

In  this  paper,  however,  the  authors  are  interested  in  obtaining 
the  entire  parameter  probability  distribution.  This  distribution 
can  be  used  to  extract  both  the  parameter  estimate  and  credible 
intervals  which  provide  a  measure  of  confidence  in  the  estimate. 
The  authors  view  well-defined  credible  intervals  as  essential  in 
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damage  detection  application.  Optimal  decisions  regarding  the 
maintenance  and  usage  of  a  structure  are  predicated  on 
quantifying  the  uncertainty  in  the  damage  estimate.  To  this  end 
the  Bayesian  estimation  framework  is  conceptually  the  most 
appealing.  In  Bayesian  estimation  one  combines  prior  information 
with  the  likelihood  to  produce  the  desired  posterior  parameter 
estimate  A  number  of  researchers  have  had  success  using  a  Mar- 
kov-Chain-Monte-Carlo  (MCMC)  approach  to  explore  the  posterior 
distribution  in  structural  parameter  estimation  problems.  Numer¬ 
ous  works  by  Beck  et  al.  (see  e.g.  [9-11])  have  used  MCMC  methods 
to  estimate  structural  parameters.  Additional  work  by  Glaser  et  al. 
[12]  illustrated  the  approach  in  detecting  stiffness  reduction  in 
beams  using  static  measurements.  Zhang  and  Cho  [13]  also  used 
the  MCMC  approach  to  help  design  an  evolutionary  algorithm  for 
performing  system  identification. 

The  MCMC  approach  (to  be  described)  works  by  repeatedly  per¬ 
turbing  the  parameter(s)  of  interest  and  either  accepting  or  reject¬ 
ing  the  proposals  based  on  a  well-defined  criteria.  However,  the 
MCMC  algorithm  can  become  “stuck”  in  a  local  minima,  again 
due  to  the  complicated  likelihood  function  (see  Section  2).  This 
is,  of  course,  why  so  many  researchers  use  GAs  to  explore  compli¬ 
cated  parameter  spaces  in  obtaining  MLEs.  As  a  result,  recent 
works  have  proposed  to  fuse  GAs  with  MCMC  in  order  to  more  effi¬ 
ciently  explore  the  posterior  distribution.  Efforts  to  this  end  in¬ 
clude  the  work  of  Vrugt  et  al.  [14]  who  combined  differential 
evolution  (a  genetic  algorithm  search  technique)  with  MCMC  to 
draw  samples  from  parameters  in  a  complex  soil  moisture  model. 
Also  noted  is  the  work  of  Zhang  and  Cho  [13]  who  combined  the 
searching  capabilities  of  an  evolutionary  algorithm  with  MCMC 
to  produce  faster  convergence  of  parameter  estimates  using  data 
generated  from  a  laser  system.  The  use  of  multiple  solutions  that 
can  exchange  information  as  they  explore  the  posterior  have  come 
to  be  known  as  “population-based  MCMC”  methods.  A  nice  over¬ 
view  of  population-based  MCMC  (Pop-MCMC)  methods  is  provided 
byjasra  et  al.  [15]  where  the  initial  idea  was  credited  to  Geyer  [16]. 
Other  recent  Pop-MCMC  works  include  the  imaging  application  of 
Kim  et  al.  [17].  Although  this  approach  has  seen  recent  use  in  these 
and  other  complicated  estimation  problems  it  has  not  yet  been 
used  in  structural  dynamics  problems  despite  obvious  advantages 
over  conventional  MCMC  in  sampling  complicated  posterior 
distributions. 

Our  goal  in  this  manuscript  is  therefore  to  adopt  the  Pop-MCMC 
approach  to  a  difficult  parameter  estimation  problem  in  structural 
dynamics,  namely  the  estimation  of  parameters  governing  crack 
damage  in  an  aluminum  plate.  The  specific  implementation  of 
the  pop-MCMC  algorithm  is  new  as  is  the  application  to  structural 
vibration  problems.  The  population-based  approach  is  shown  to 
provide  reliable  estimates  of  the  crack  parameters  under  realistic 
levels  of  noise.  These  estimates  are  not  always  possible  using  stan¬ 
dard  MCMC  as  will  be  shown.  Additionally,  implementation  of  the 
Pop-MCMC  algorithm  requires  that  the  forward  model  be  compu¬ 
tationally  efficient.  Section  4  is  therefore  devoted  to  developing  a 
low-dimensional  finite  element  (FE)  crack  model.  The  model  uses 
specially  tailored  elements,  resulting  in  a  model  that  retains  the 
fidelity  of  the  standard  FE  approach  without  the  computational 
overhead. 


s  =  w(0)+£  (1) 

where  the  noise  vector  {={,■,  i=  1 . T  is  comprised  of  T,  zero- 

mean,  iid  entries  with  joint  probability  distribution 


Ph(0 


l  -E^2 

(27l(72)r/2 


(2) 


Thus,  the  noise  is  assumed  to  be  stationary,  Gaussian  with  variance 
a2.  The  likelihood  function  is  then  formed  by  simply  substituting 
it  =  Si  -  Wj(0)  giving 


Ph(S)  s  P„(s|0) 


T 

l  -E  (si-wi(8))2/2ff2 

(2na2)T/2e 
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which  is  the  probability  of  having  observed  the  sequence  of  data  s 
given  the  model,  defined  by  parameters  0.  The  aforementioned 
MLEs  are  obtained  by  maximizing  Eq.  (16).  However,  the  goal  is 
to  obtain  estimates  of  both  the  parameter  and  the  amount  of  uncer¬ 
tainty  in  the  parameter  estimate.  For  this  reason  the  Bayesian  esti¬ 
mation  philosophy  is  adopted  thus  providing  estimates  of  the  entire 
parameter  posterior  distribution  from  which  credible  intervals  for 
our  parameters  can  easily  be  obtained. 

Bayes’  rule  states  that  the  joint  posterior  distribution  of  our 
model  parameters  may  be  found  by  the  relationship 


Pe(0) 


PH(s|0)Pn(0) 

Pd(s) 


(4) 


where  pn(0)  is  the  joint  prior  parameter  distribution,  reflecting  any 
a  priori  information  one  might  have  about  our  parameter,  and  pD( s) 
is  the  joint  distribution  of  our  acquired  data.  Of  course  the  above 
expression  is  for  the  joint  parameter  distribution  whereas  one  is 
typically  interested  in  the  marginal  (individual)  posteriors,  i.e. 
P0p(dp).  Analytically  this  would  require  integrating  Eq.  (4)  over 
the  other  P  -  1  parameters,  e.g. 

Ve p{SP)=  [  Pe(0)d0-p  (5) 

J«p-' 

where  the  notation  d0_p  denotes  the  multi-dimensional  inte¬ 
gral  over  all  parameters  other  than  0P.  This  cannot  typically  be  done 
in  closed  form  and  instead  the  authors  resort  to  numerical  methods. 

The  Markov  Chain  Monte  Carlo  (MCMC)  algorithm  was  pro¬ 
posed  by  Hastings  [18]  as  a  means  of  drawing  samples  from  Eq. 
(4)  directly  (without  having  to  integrate).  The  MCMC  algorithm  is 
now  fairly  standard  and  the  reader  is  referred  to  [19]  for  a  descrip¬ 
tion  in  the  context  of  structural  dynamics.  The  algorithm  builds  a 
Markov  chain  for  each  parameter  such  that  samples  from  the  chain 
are,  in  fact,  samples  from  the  desired  posterior  distribution.  The 
chain  is  formed  by  first  proposing  state  transitions  for  each  param¬ 
eter  via  a  so-called  proposal  distribution,  q(Q'\0p).  This  distribution 
provides  a  rule  for  generating  a  candidate  parameter  value  0*  given 
the  current  value  in  the  chain  0P.  The  proposed  parameter  (keeping 
all  other  parameters  fixed  at  their  current  value)  is  then  accepted 
or  rejected  with  probability 


rMH  =  min 


ph(s\q;)Pup{s;)  l 
PH(S|0p)Pnp(0p)’  J 


(6) 


2.  Brief  review  of  Bayes  rule  and  MCMC 

Given  data  from  a  system  of  interest,  denoted  by  the  vector  of  N 
observations  s,  and  a  model  that  describes  those  data  w (0),  the 
practitioner’s  job  is  to  estimate  the  model  parameter  vector 

0  =  (0lt02 . Op).  Typically  one  takes  “good"  estimates  to  be  those 

that  perform  well  in  the  face  of  the  inevitable  uncertainty  (e.g. 
noise)  in  the  data.  In  this  paper  it  is  assumed  that  the  uncertainty 
takes  the  form  of  additive  noise  such  that  our  measured  data  are 


This  procedure  is  referred  to  as  the  Metropolis-Hastings  (MH)  algo¬ 
rithm  after  the  originators  [20,18].  The  same  is  done  for  each  of  the 
p  =  1,  . .  .,P  parameters  in  0  while  holding  the  other  P  -  1  parameters 
fixed.  Thus,  the  approach  is  really  drawing  samples  from  the  condi¬ 
tional  posterior  Pgp(0p |0_p).  In  other  words,  given  that  the  other 
parameters  are  fixed  to  their  current  values  in  the  chain,  sample 
the  posterior  for  parameter  p.  This  procedure  is  referred  to  as  Gibbs 
sampling  and  eliminates  the  need  to  perform  the  high-dimensional 
integral  required  by  Eq.  (5).  The  above-described  procedure 
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constitutes  one  iteration  of  the  Markov  chain.  Starting  with  initial 
parameter  values  chosen  from  their  respective  prior  distributions 
0P( 0)  ~  Pnr(8p)  this  process  (accept/reject)  repeats  I<  times  yielding 
the  samples  0p(k),  k  =  1 ...  ,/C 

The  distribution  q(0*|0p)  can  be  thought  of  as  a  rule  for  perturb¬ 
ing  0P.  The  magnitude  of  this  perturbation  will  be  decided  by  the 
parameters  associated  with  q().  For  example,  it  is  common  to 
choose  the  proposal  distribution  to  be  the  Uniform  distribution 
centered  at  0P,  i.e. 

q(e;\ep)  =  U(0p~A.8p+A)  (7) 

where  A  is  a  real  constant  that  plays  the  role  of  the  perturbation 
size.  Typically  A  is  tuned  dynamically  so  that  on  average 
rMH  =  03  -  0.5  (see  [19]  for  a  more  thorough  discussion  and  for 
sample  code).  Heuristically  it  is  easy  to  see  how  the  approach 
works.  The  algorithm  continually  perturbs  each  parameter,  checks 
whether  or  not  a  better  fit  (consistent  with  the  prior)  is  achieved 
by  computing  rMn,  and  keeps  parameter  values  that  do  well  in  this 
regard.  The  fact  that  in  the  end  the  chain  of  values  are  samples  from 
Pgp(0p)  is  actually  quite  remarkable.  However,  one  can  immediately 
see  where  problems  can  arise. 

Consider  a  bimodal  posterior  distribution  of  a  single  parameter 

0, 


is  the  true  posterior  distribution.  In  what  follows  the  authors  use 
n  =  1  to  denote  the  true  posterior  distribution  from  which  the  de¬ 
sired  samples  will  be  drawn.  If  an  irreducible,  a-periodic  Markov 
chain  can  be  constructed  that  has  pc(0(JNi)  has  its  invariant  distri¬ 
bution,  then  samples  from  the  marginal  distribution 

Pe(fl)  =  /  Pc(0(,:N))d0(2:N) 

JuN-' 

can  be  drawn  where  the  notation  JrN- ,  d0(2:N)  denotes  the  multi¬ 
dimensional  integral  over  all  parameters  vectors  other  than 
0  =  0(1).  This  is  accomplished  numerically  by  running  “N"  Markov 
chains  concurrently,  each  exploring  its  own  posterior  distribution 
p„(0(n)).  Each  chain  can  be  considered  in  turn,  holding  the  others 
fixed  (Gibbs  sampling),  and  samples  from  p i(0(1))  are  retained  as 
samples  from  the  desired  joint  posterior. 

There  exists  a  good  deal  of  freedom  in  choosing  the  pn(  )  and  in 
selecting  different  types  of  proposals  for  exploring  the  parameter 
space.  With  regard  to  the  former  one  would  like  to  wisely  choose 
the  p „(■)  so  as  to  facilitate  easy  exploration  of  the  parameter  space. 
While  other  approaches  can  be  used  (see  [15]),  in  this  paper  the  so- 
called  tempered  sequence  of  distributions  is  used 

pn(0<">)=pg(0‘">)  (10) 


p^=^e^2/(2fa2) 

I  1  „-(e+rt2/(2(i-tV)  (g) 

2^271(1  -f)a2 

for  known,  positive  constants/,  a,  p.  The  constant/<  1  specifies  the 
fraction  of  the  distribution  variance  associated  with  the  distribution 
peaks  at  +p,  -p.  For  this  example  the  values  a  =  l,/=  0.3,  were  used 
along  with  two  different  values  for  the  peak  locations,  p  =  2,4.  The 
resulting  distribution  presents  difficulties  for  the  conventional 
MCMC  algorithm  for  the  p  =  4  case  as  shown  in  Fig.  1.  For  p  =  2 
the  peaks  of  the  distributions  are  close  enough  that  the  algorithm 
can  easily  move  back  and  forth  between  the  two  regions  of  high 
probability.  However  for  p  =  4  the  Markov  Chain  quickly  becomes 
“trapped”  in  a  single  portion  of  the  posterior  distribution.  This  is 
simply  due  to  the  fact  that  the  proposal  distribution  is  not  capable 
of  moving  the  chain  easily  from  one  peak  to  the  other.  The  solution, 
however,  is  not  as  simple  as  changing  the  proposal  distribution  to 
give  us  larger  perturbations  to  our  chain  (i.e.  increase  A).  If  the  pro¬ 
posal  values  are  very  far  from  the  existing  values  they  will  almost 
always  be  rejected,  thus  the  Markov  Chains  will  take  an  extremely 
long  time  to  converge.  An  efficient  sampler  is  one  that  would  allow 
us  to  locally  explore  a  high  probability  region  of  the  posterior  while 
simultaneously  provide  a  mechanism  for  covering  large  distances  in 
parameter  space  to  reach  other  high  probability  regions.  The  popu¬ 
lation-based  MCMC  approach,  described  next,  was  designed  specif¬ 
ically  for  this  reason. 

3.  Population-based  MCMC 

As  with  standard  MCMC,  the  goal  of  the  Pop-MCMC  algorithm  is 
to  draw  samples  from  some  desired  posterior  distribution  p®(0).  A 
very  nice  introduction  to  the  topic  of  Population-based  MCMC 
(Pop-MCMC)  is  given  by  Jasra  et  al.  [15].  The  basic  idea  is  to  first 
create  a  new,  composite  posterior  density 

pc(0(1:N))  =  f[pn(0(n))  (9) 

n-l 

which  is  a  function  of  the  composite  parameter  vector 
0(|:N)  =  (0(1), . . .  ,0(N)).  It  is  required  that  p„(0(n))  =  p@(0)  for  at  least 
one  n,  i.e.  one  of  the  posteriors  comprising  the  composite  density 


for  ( n  e  (0, 1  ].  For  (h  =  1,  of  course,  one  has  the  true  posterior.  For 
smaller  values  of  £„  one  obtain  successively  smoother  versions  of 
the  original  posterior.  The  idea  is  that  the  smoothed  distributions 
are  easier  for  their  corresponding  Markov  chains  to  explore,  yet 
they  are  still  related  to  the  true  posterior  and  therefore  still  carry 
information  about  high  probability  regions  of  the  parameter  space. 
As  an  example  consider  N  =  4  separate  instances  of  the  bi-modal 
distribution  (8)  raised  to  the  £ „  =  1,  0.5,  0.1,  0.0125  powers,  respec¬ 
tively.  These  distributions  are  shown  in  Fig.  2.  The  standard  MCMC 
algorithm  will  have  an  easier  time  exploring  these  smoother 
parameter  spaces.  Of  course  only  samples  from  the  distribution 
with  Ci  =  1  are  of  interest,  however  the  additional  distributions 
can  clearly  facilitate  efficient  sampling  if  one  has  a  means  of  passing 
information  between  chains.  The  chain  associated  with  the  true 
posterior  needs  to  be  informed  by  the  chains  exploring  the  smooth¬ 
er  distributions. 

Perhaps  the  simplest  type  of  move  to  accomplish  this  is  the  so- 
called  swap  move.  This  is  similar  to  a  Metropolis-Hastings  move 
where  the  proposal  is  to  consider  swapping  the  parameter  values 
in  chains  u,  v.  Assuming  an  equal  probability  of  selecting  chains 
u,  v  from  the  N  possibilities,  this  swap  is  accepted  with  probability 


t swap 


min 


Pu(9(>,(ew)  A 
Pu(o(u))pAom)’  I 


(ii) 


Typically  swap  moves  are  not  proposed  after  every  iteration  in  the 
Markov  Chains  but  are  performed  with  some  probability.  Returning 
to  the  bi-modal  example,  consider  N  =  4  chains  exploring  the  com¬ 
posite  target 

Pc(0(,:4))=lW',)).  (12) 

n-l 


thus  one  has  four  Markov  chains  running  concurrently  with  £„  =  1, 
0.75,  0.5,  0.25,  respectively.  After  each  iteration  in  the  Markov  chain 
a  swap  move  is  performed  with  50%  probability.  This  move  consists 
of  uniformly  selecting  two  of  the  chains  and  evaluating  Eq.  (11).  If 
the  proposal  is  accepted  the  values  in  the  chains  are  exchanged 
and  the  algorithm  continues  to  the  next  iteration  in  the  Markov 
chains.  Fig.  3  shows  the  results  of  this  sampler  All  four  Markov 
chains  are  informing  each  other  as  to  the  presence  of  multiple  peaks 
in  the  distribution.  The  end  result  is  that  the  chain  associated  with 
£ i  =  l  contains  samples  from  the  desired  posterior  distribution;  this 
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Fig.  1.  Markov  Chains  and  estimated  bi-modal  posterior  distribution  given  by  Eq.  (8)  with  \x  =  2  (a,b)  and  [i  =  4  (c,d),  respectively. 


Fig.  2.  Successive  p„(0)  =  p@(0)  corresponding  to  the  bi-modal  posterior  distribu¬ 
tion  Eq.  (8).  The  true  posterior  distribution  from  which  samples  are  sought  is  given 

by  c,  =  i- 


distribution  is  shown  in  Fig.  3b  and  compares  favorably  to  the  true 
posterior  density. 

While  the  swap  move  is  an  effective  means  of  communicating 
among  chains,  a  more  sophisticated  type  of  move  is  required  for 
the  structural  estimation  problem  considered  next.  Particularly  for 
multivariate  parameter  estimation,  it  is  useful  to  design  a  move  that 


can  hold  multiple  parameters  fixed  while  allowing  other  groups  of 
parameters  to  move  simultaneously.  Such  a  move  is  not  allowed 
with  the  standard  Metropolis-within-Gibbs  sampling  (one  parame¬ 
ter  moved  at  a  time)  or  by  the  swap  move  which  moves  all  parame¬ 
ters  simultaneously.  A  particularly  effective  move  is  to  use 
differential  evolution  to  generate  a  trial  vector.  Differential  evolution 
(DE)  is  the  engine  of  a  popular  GA  used  in  searching  complex  param¬ 
eter  space  [21  [.The  approach  draws  at  random  three  members  of  the 
population,  0(u),  0(v\  0(w)  and  generates  the  trial  vector 

0(u'>  =  0<")  +  y(0W  -  gW)  (13) 

where  y  is  a  user-defined  constant.  Each  of  the  P  elements  in  this 
trial  vector  replace  the  elements  of  the  original  vector  0(u)  with 
50%  probability.  This  final  step  (keep  new  element  or  retain  old) 
emulates  the  “cross-over”  step  common  to  most  GAs.  Once  the  trial 
vector  has  been  generated  it  is  accepted/rejected  using  rMH  (Eq.  6). 
However,  in  order  to  differentiate  among  the  types  of  moves  the 
acceptance  ratio  for  the  DE  move  will  be  denoted  rDE. 

It  should  mentioned  that  others  have  proposed  a  straight  cross¬ 
over  move  whereby  the  parameter  vector  is  split  at  a  point  m  <  P 
for  two  “parent”  vectors  and  generating  the  trial  vectors: 

It*  =  ..,<#>] 

and  accepting  with  probability  r^p  (Eq.  1 1 )[  1 5 ].  This  is  indeed  a 
useful  way  of  exchanging  information  between  chains,  however 
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Fig.  3.  (a)  Markov  chains  and  (b)  estimated  posterior  distribution  associated  with  the  bi-modal  distribution  given  by  Eq.  (8). 


the  authors  have  found  that  this  mechanism  is  already  provided  for 
in  the  differential  evolution  move,  thus  this  type  of  move  is  not  used 
in  this  paper. 

The  above-described  moves  clearly  borrows  from  the  genetic 
algorithm  (GA)  approach  to  optimization  problems.  In  fact,  Pop- 
MCMC  is  effectively  combining  the  efficient  search  capabilities  of 
GAs  with  the  power  of  the  Bayesian  MCMC  approach  to  sampling. 
Both  types  of  moves,  swap  and  differential  evolution,  will  be  used 
in  the  structural  dynamics  example  presented  in  subsequent 
sections. 


4.  An  efficient  cracked  plate  model 

One  critical  component  of  this  process  is  a  computationally  fast 
model  that  is  also  flexible  enough  to  incorporate  a  variety  of  dam¬ 
age  modes.  The  model  predicts  the  lateral  deflection  of  the  plate  at 
a  given  point  in  space  and  time,  w(x,y,  t).  Efficiency  of  the  forward 
model  is  essential  as  the  Pop-MCMC  algorithm  requires  a  very 
large  number  of  model-to-date  comparisons,  i.e.  evaluate  rMH.  rswap, 
rDE.  Flexibility  is  also  desired  because  while  the  current  study  con¬ 
siders  damage  consisting  of  only  a  single  straight  crack,  it  is  the 
model  alone  that  limits  the  type  of  damage  that  can  be  considered 
by  this  technique.  A  model  that  could  handle  multiple  cracks, 
branched  cracks,  edge  cracks,  or  corrosion  would  be  valuable.  A 
purely  analytical  model  entailing  virtually  no  computation  time 
would  be  ideal  on  both  counts.  Unfortunately,  analytical  solutions 
for  the  dynamic  response  of  a  rectangular  plate  with  an  arbitrary 
crack  (position,  orientation,  and  size)  do  not  presently  exist  to 
the  knowledge  of  the  authors.  Solecki  [22]  produced  a  solution 
for  the  natural  frequencies  of  a  plate  with  a  single  crack,  but  did 
not  extend  it  to  the  myriad  of  other  situations  of  future  interest. 
A  second  possible  family  of  models  are  finite  element  solutions. 
These  clearly  have  the  desired  flexibility,  but  involve  time  march¬ 
ing  and  so  are  computationally  intense.  In  this  present  work,  a 
compromise  approach  is  taken.  A  finite  element  eigen-solution  is 
used  to  build  in  the  features  of  the  crack  singularity.  Building  off 
of  these  numerical  frequencies  and  mode  shapes,  modal  analysis 


is  used  to  describe  the  actual  motion  in  space  and  time,  w(x,y,t). 
The  eigen-solution  to  this  finite  element  model  gives  the  modes 
in  terms  of  the  nodal  displacements.  The  values  at  the  instrumen¬ 
tation  locations  are  interpolated  via  the  natural  neighbor  method. 
Standard  modal  analysis  proceeds  from  there.  This  hybrid  solution 
is  many  times  faster  than  time  marching,  while  retaining  the  flex¬ 
ibility  inherent  in  the  finite  element  approach. 

An  in-house  finite  element  code  was  developed  for  this  partic¬ 
ular  application,  so  as  to  more  easily  wrap  the  MCMC  process 
(code)  around  the  time  series  solution.  The  model  parameters  are 
the  location  of  the  center  of  the  crack  (xcrack,ycrack),  the  crack  length 
(a),  and  the  orientation  of  the  crack  measured  from  the  positive  x- 
axis  (a).  See  Fig.  4c.  The  elements  away  from  the  crack  tip  are  stan¬ 
dard  eight-noded  quadrilateral  Mindlin  serendipity  elements,  with 
nodes  at  each  corner  and  at  the  middle  of  each  side.  Adjacent  to  the 
crack  tip  the  eight-noded  quads  are  modified  as  described  in  Refs. 
[23,24].  In  short,  two  corner  nodes  and  the  node  in  between  them 
are  moved  to  a  single  location,  leaving  a  triangular  element.  The 
adjacent  side  nodes  are  then  moved  from  the  midpoints  to  one 
quarter  of  the  side  length.  The  result  is  an  element  with  a  stress 
field  that  varies  as  Jp  where  r  is  the  distance  from  the  collapsed 
node  (to  be  placed  on  the  crack  tip).  See  Fig.  4a  and  b  for  an  illus¬ 
tration.  This  stress  distribution  exactly  matches  the  Mitchell  solu- 
tion[25],  which  is  the  analytic  stress  field  near  a  static  crack  tip. 
Because  these  augmented  triangular  elements  capture  the  crack 
tip  behavior,  they  are  placed  around  the  crack  tip  in  a  pinwheel 
fashion.  See  Fig.  4d.  This  permits  a  more  sparse  mesh  near  the 
crack  tip  than  would  otherwise  be  necessary.  This  reduces  the 
number  of  degrees  of  freedom  required  and,  hence,  the  computa¬ 
tion  time. 

The  cracked  plate  model  involves  several  standard  assump¬ 
tions.  First,  the  model  is  linear  and  it  is  assumed  that  material 
properties  are  known  exactly  (though  these  could  be  left  as  un¬ 
known  parameters  and  found  via  the  Bayesian/MCMC  approach). 
The  deflections  are  presumed  small  and  crack  growth  is  not  con¬ 
sidered;  the  latter  is  reasonable  under  small  deflections.  The  crack 
is  also  assumed  to  remain  open,  such  that  impacts  at  the  crack 
interface  are  ignored.  Any  mass  lost  loss  due  to  the  crack  is 
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Fig.  4.  (a)  Nodal  configuration  of  the  Mindlin  serendipity  element,  (b)  resulting  interpolation  functions,  (c)  schematic  of  plate  with  crack  parameters  and  (d)  resulting  course 
FEM  mesh. 


presumed  negligible.  Finally,  there  is  only  one  crack  and  the  tips 
are  at  least  one  half  of  the  crack  length  from  the  edges  of  the 
plate.  Since  finding  smaller  cracks  is  the  more  interesting  and 
useful  problem,  the  allowable  zone  is  still  a  very  large  fraction 
of  the  plates  area.  Even  taken  together,  the  authors  do  not  believe 
that  these  requirements  are  too  onerous. 

As  with  any  finite  element  analysis,  convergence  of  the  mesh 
must  be  verified.  Due  to  the  large  number  of  iterations,  it  is  not 
practical  to  verify  convergence  for  every  perturbed  set  of  parame¬ 
ters.  Instead,  because  all  the  meshes  used  are  qualitatively  similar, 
convergence  was  checked  on  several  representative  parameter 
vectors. 


T  M 

Q(M)  =  -  w>/(0))2, 

i  j 


the  sequence  of  tempered  posterior  distributions  is  formed  as 


Pn(e^)  =  pH(s\e^fpn(e^) 


n  =  1,..., N 


(2nalf™/2  " 


(17) 


using  the  sequence  (h  =  1.0, 

‘sIJ+l  n  =  1,  — N  -  1 


5.  Implementation 

In  this  paper  the  structure  of  interest  assume  a  clamped  rectan¬ 
gular  plate  measuring  1 .25  m  x  1  m  with  thickness  h  =  0.01  m  with 
material  properties  E  =  209  GPa  (Young’s  modulus),  v  =  0.3  (Pois¬ 
son’s  ratio),  and  p  =  7850  kg/m3  (density).  It  will  also  be  assumed 
that  the  plate  has  been  instrumented  with  M  displacement  sensors 
capable  of  sampling  the  plate’s  response  to  an  input  at  /  =  1,. .  .,T 
equally  spaced  points  in  time.  The  observed  signal  model  is  there¬ 
fore  written  (as  before) 

Sij  =  w,j(0)  +  £,j,  i=\,...T,  (15) 

where  the  are  taken  as  realizations  of  an  iid,  zero-mean,  Gaussian 
random  process,  i.e.  each  qj  ~  N(0,cr2).  The  model  (described  in  the 

previous  section)  is  therefore  evaluated  at  the  same  times  i  =  1 . T 

and  locations  j=l _ ,M  as  the  observed  data.  Under  this  noise 

model  the  likelihood  function  for  the  data  is  given  by  [19] 


Ph(s|0) 


T  M 

(2ncr2)™'2e 


Denoting  the  sum-squared  error  over  all  sensors  by 


(16) 


as  suggested  byjasra  et  al.  [15].  The  idea  here,  as  in  the  toy  exam¬ 
ple,  is  to  explore  a  composite  posterior  where  the  “smoothed"  mar¬ 
ginal  posteriors  n  =  2,...,N  are  related  to  the  true  posterior  of 
interest,  pi(0(1)). 

This  implementation  of  the  pop-based  algorithm  proceeds  as 
follows.  The  parameter  vectors  for  each  of  the  N  chains  are  initial¬ 
ized  by  drawing  samples  from  the  priors.  Then,  for  each  iteration  in 
the  Markov  chain,  one  of  the  chains  n  g  N  is  selected  with  uniform 
probability,  and  a  standard  MH  update  is  performed  for  each  of  the 
P  parameters  in  0<n)  using  the  Gibbs  sampling  strategy.  Thus,  for 
each  parameter  p  =  1,. .  .,P  one  generates  a  candidate  value  using 
Eq.  (7),  evaluates  the  ratio 


rMH  =  Exp 


Cn 

2<t2 


(Q(s,0;(")|0Ln(J)~Q(s,0(">|0(">)) 


P,(Q 

p*(fiD 


and  accepts  with  probability  min(rMH.l).  The  (unknown)  noise  var¬ 
iance  associated  with  the  nth  chain,  an,  also  needs  to  be  sampled.  It 
has  been  demonstrated  [19]  that  by  choosing  a  vague  prior  for  this 
parameter,  one  may  directly  sample  from  the  posterior  via 


r(MT/2, 2/(£„Q.(s,  0(n>)) 


where  F(a,b)  denotes  the  Gamma  distribution  with  parameters  a,b. 
Once  each  of  the  parameters  (including  the  noise  variance)  have  been 
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updated  for  chain  n,  two  different  chains  u,  ve  [1,N]  are  selected  at 
random  with  uniform  probability.  With  50%  probability  either  a 
“swap”  move  or  a  “DE”  move  between  these  chains  is  performed. 
For  the  swap  move,  the  prior  ratios  will  cancel,  thus  one  evaluates 


=  Exp 


— 4u 


jhQM-’l-jhQM") 


(IS) 


accepting  the  move  with  probability  min(rswap,  1).  For  the  DE  move 
one  requires  three  randomly  drawn  chains  (see  Eq.  13)  to  generate 


the  trial  vector  0,u'\  The  ratio  of  priors  in  the  acceptance  criteria  is 
also  required  as  they  do  not  cancel  out.  One  therefore  evaluates 


rDE  =  Exp 


^(q(s,0(uVQ(M(u))) 


n  fpPA^]) 
n  pP«(C)  ’ 


(19) 


and  accept  the  move  with  probability  min(rD£,  1 ).  To  summarize,  the 
algorithm  picks  a  chain  at  random  and  performs  a  standard  MH  up¬ 
date  on  each  of  the  parameters,  including  the  noise  variance.  The 
algorithm  then  selects  chains  at  random  and  performs  either  a  swap 
move  (Eq.  18)  or  crossover  (Eq.  19)  with  50%  probability.  This 


Fig.  5.  Main  argument  of  the  likelihood  function,  Q( s,6),  plotted  as  a  function  of  (a)  crack  location  on  the  plate  and  (b)  the  parameters  a  (crack  length)  and  6C,  crack 
orientation.  For  each  plot  the  remaining  parameters  were  held  fixed  at  their  true  values.  True  parameter  values  are  denoted  with  an  “X”. 
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Fig.  6.  Estimated  posterior  densities  associated  with  (a)  crack  length,  (b,c)  crack  location,  and  (d)  crack  orientation.  The  true  parameter  values  are  a  =  0.1,  (xc,yc)  =  (0.3, 0.5), 
and  0C  =  30. 
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procedure  repeats  for  some  number  of  iterations  until  enough  sam¬ 
ples  have  been  drawn  from  the  posterior  distribution. 

The  algorithm  is  fairly  simple  to  implement  in  software,  how¬ 
ever  is  clearly  computationally  intensive.  Each  evaluation  of  a  ratio 
(cMwCswap  or  cDe)  requires  solving  the  forward  model  described  in 
Section  4.  The  efficiency  of  the  forward  model  is  therefore  of  great 
importance  in  using  MCMC  in  structural  system  identification 
problems. 


6.  Cracked  plate  identification 

There  are  four  parameters  that  determine  the  state  of  damage  in 
the  plate  model:  crack  length  a,  crack  location  (xc,yc),  and  crack  ori¬ 
entation  0C.  These  damage  parameters  are  fixed  to  the  values 
a  =  0.1  m,  xc  =  0.3  m,  yc  =  0.5  m,  and  8C  =  30°.  Further  assume  that 
four  displacement  sensors  have  been  placed  on  the  surface  of  the 
plate  at  the  x-y  locations  (0.375,0.375),  (0.375,0.862),  (0.862,0.37 
5),  (0.862,0.862).  The  acquired  data  will  consist  of  each  sensor’s  re¬ 
sponse  to  four  separate  impacts  (hammer  strikes)  at  locations 
(0.29, 0.275),  (0.29, 0.725),  (0.96, 0.275),  and  (0.96, 0.275).  These  sen¬ 
sor  locations  were  chosen  so  as  to  maximize  the  sum  of  the  first  four 
modes  in  the  response.  No  claims  of  optimality  are  made  regarding 
this  choice,  in  fact  finding  sensor  locations  that  produce  a  more  well- 
defined  likelihood  is  an  active  area  of  research. 

The  data  used  in  the  identification  procedure  consist  of  four 
simulated  impulse  response  signals,  sampled  at  4  kHz  for  a 


duration  of  2  s  (T  =  8000  observations).  The  first  nine  modes  were 
used  in  generating  the  solution,  all  well  below  the  Nyquist  fre¬ 
quency  of  2  kHz.  The  signal  to  noise  ratio  was  set  at  17  dB  (50:1) 
corresponding  to  the  level  of  noise  observed  in  previous 
experiments. 

Given  this  response  the  goal  is  to  estimate  the  posterior  distri¬ 
butions  associated  with  each  of  the  crack  parameters.  Using  the 
standard  MCMC  algorithm  the  authors  were  unable  to  consistently 
identify  the  parameters  due  to  the  aforementioned  problem  with 
exploring  complicated  likelihood  functions.  Fig.  5  shows  the  argu¬ 
ment  of  the  likelihood,  Q(s,0 )  as  a  function  of  different  parameter 
combinations  (holding  the  others  fixed  at  their  true  values)  As  with 
the  simple  1-D  example  (Fig.  1),  one  can  immediately  see  the  dif¬ 
ficulty  presented  by  this  estimation  problem.  Consider  Fig.  5a. 
The  location  parameter  xc  has  a  minimum  at  the  true  location, 
xc  =  0.3,  however,  there  is  also  a  “trough”  of  minima  with  a  partic¬ 
ularly  low  value  at  (xc  =  0.85,  yc  =  0.28).  A  similarly  complicated 
likelihood  is  shown  for  the  parameters  a,  9C  (Fig.  5b).  These  plots 
illustrate  the  fundamental  difficulty  in  parameter  identification 
for  structural  systems:  multiple  damage  states  can  yield  very  similar 
structural  response  data. 

For  this  reason  the  authors  have  gravitated  toward  the  popula¬ 
tion-based  approach  to  sampling  the  posterior.  In  the  example  that 
follows,  N  =  10  chains  were  used  for  each  of  the  parameters  a,  0C,  xc, 
yc.  The  prior  for  the  crack  parameter  a  is  take  to  be  a  Gamma  dis¬ 
tribution,  pIa(a)  =  a*-1  with  parameters  a.  =  1.25,  p  =  0.025. 

Thus,  essentially  no  damage  is  assumed  at  the  outset.  Certainly  this 
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Fig.  7.  Estimated  posterior  densities  associated  with  (a)  crack  length,  (b,c)  crack  location,  and  (d)  crack  orientation.  The  true  parameter  values  are  a  =  0.1,  (xc,yc)  =  (0.6,0.35), 
and  9C  =  -20. 
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distribution  can  be  altered  to  reflect  a  known  level  of  damage, 
however  the  more  typical  case  is  to  assume  the  plate  is  healthy 
and  allow  the  data  (likelihood)  to  drive  the  posterior.  The  form  of 
this  prior  (Gamma)  was  chosen  based  on  the  fact  that  the  crack 
cannot  be  negative.  One  might  also  have  chosen  a  Beta  prior  which 
has  finite  support  on  the  plate,  however  for  the  prior  parameters 
chosen  there  is  essentially  no  probability  of  the  crack  extending 
off  the  end  of  the  plate.  For  the  crack  location  parameters  xc,  yc  uni¬ 
form  priors  were  chosen  over  the  span  of  the  crack,  i.e.  xc  U( 0, 1 .25 ), 
yc  1/(0, 1 )  reflecting  the  fact  that  no  a  prior  knowledge  about  the 
crack  location  exists.  Similarly  the  prior  on  crack  angle  was  taken 
as  0C  U(-90,90)  as  in  general  one  will  not  know  the  crack  orienta¬ 
tion  a  priori  either.  The  prior  distributions  are  displayed  in  Fig.  6. 

Initializing  the  Markov  chains  using  these  priors,  the  popula¬ 
tion-based  MCMC  algorithm  was  run  for  60,000  iterations  with  a 
burn-in  of  50,000  iterations.  The  remaining  10,000  iterations  were 
stored  and  used  to  form  the  posterior  densities  shown  in  Fig.  6.  In 
this  simulation,  and  in  others  the  authors  have  looked  at,  crack 
length  is  perhaps  the  most  easily  identified  parameter.  The  Markov 
chains  in  all  populations  tend  to  the  true  value  (a  =  0.1 )  after  only  a 
few  thousand  iterations.  Likewise,  the  location  parameter  yc  exhib¬ 
its  a  uni-modal  posterior  distribution.  The  other  parameters,  xc,  dc 
are  significantly  more  challenging  to  estimate  as  both  show  the 
clear  presence  of  multiple  maxima  in  the  likelihood.  The  location 
parameter  xc  exhibits  multi-modal  behavior.  In  particular  one  can 
see  the  multiple  minima  observed  in  Fig.  5.  The  two  peaks  located 
at  xc  =  0.3,  xc  =  0.4  presented  a  significant  challenge  for  the  regular 
MCMC  algorithm.  Because  crack  lengths  were  initially  assumed 
small,  the  first  minimum  encountered  by  the  Markov  chain  was 
the  one  located  at  xc  =  0.4  and  the  solution  would  often  remain 
here.  Running  multiple  chains  easily  overcomes  this  problem.  Per¬ 
haps  the  greatest  utility  of  the  approach,  however,  can  be  seen  in 
the  identified  orientation  0C.  Regardless  of  crack  geometry,  there 
are  nearly  always  multiple  solutions  for  0C  that  come  close  to  max¬ 
imizing  the  likelihood.  The  result  are  multiple  well-defined  local 
maxima  that  easily  trap  the  standard  MCMC  algorithm.  By  con¬ 
trast,  Fig.  6  captures  the  relative  heights  of  these  maxima  indicat¬ 
ing  that  the  highest  probability  for  crack  orientation  is  at  the  true 
value,  0C  =  30°. 

As  a  second  example,  Fig.  7  shows  the  identified  posterior  dis¬ 
tributions  for  the  case  where  the  true  crack  parameters  were  set 
to  the  values  a  =  0.1,  xc  =  0.6,  yc  =  0.35,  0C=  -20.  Additionally,  the 
Gamma  prior  (biased  toward  no  crack)  was  changed  to  a  Uniform 
prior  in  order  to  demonstrate  the  insensitivity  of  the  approach  to 
this  choice.  It  could  be  that  the  practitioner  has  no  a  prior  informa¬ 
tion  regarding  the  presence  and  length  of  a  crack,  thus  a  uniform 
prior  would  be  appropriate.  All  parameters  are  again  correctly 
identified,  provided  that  the  final  estimate  is  taken  as  the  maxi¬ 
mum  a  posterior  value.  Again,  one  sees  multiple  “good”  solutions 
which  can  often  trap  the  standard  MCMC  algorithm.  For  example, 
the  crack  location  parameter  yc  has  a  fairly  high  probability  of 
being  yc  =  0.7  despite  the  fact  that  the  true  value  is  a  factor  of 
two  different.  Similarly  multiple  peaks  for  the  orientation  parame¬ 
ter  0C  can  be  seen.  The  authors  note  that  the  standard  MCMC  algo¬ 
rithm  was  not  able  to  sample  from  these  complex  posterior 
distributions. 

The  complexity  of  the  posterior  distributions  is  simply  a  con¬ 
sequence  of  trying  to  identify  parameters  that  have  little  affect 
on  the  global  vibrations.  It  is  obvious  from  the  likelihood  plots 
(Fig.  5)  that  varying  crack  configurations  can  lead  to  nearly  the 
same  vibrational  response.  This  is  part  of  the  physics  of  structural 
identification  problems  and  the  Bayesian  approach,  using  the 
pop-MCMC  algorithm,  is  correctly  reflecting  this  uncertainty. 
Rather  than  producing  point  estimates  which  may  or  may  not 
coincide  with  the  true  parameter  value,  the  approach  described 
here  allows  for  a  full  accounting  of  the  probability  associated 


9 

with  different  crack  configurations.  The  authors  view  this  as  a 
strength  in  difficult  system  identification  problems,  particularly 
those  in  damage  detection  where  quantifying  uncertainty  is  crit¬ 
ical  to  making  informed  decisions  regarding  the  health  of  a  struc¬ 
ture.  Choosing  to  use  more  modes  in  the  model  helps  the 
situation,  however,  in  practice  it  becomes  increasingly  difficult 
to  match  model  to  data  for  higher  modes.  In  this  paper  the  focus 
was  on  the  first  five  modes  of  the  response. 

7.  Conclusions 

This  paper  proposes  a  population-based  Markov  Chain  Monte 
Carlo  method  for  solving  complicated  system  identification  prob¬ 
lems  in  structural  dynamics.  Both  the  application  and  specific 
implementation  of  the  method  are  new.  The  Pop-MCMC  method 
is  well-suited  to  the  types  of  posterior  parameter  distributions  of¬ 
ten  found  in  structural  dynamics  problems.  These  distributions 
tend  to  be  multi-modal  and  can  easily  confuse  standard  MCMC  ap¬ 
proaches.  The  population-based  approach,  by  contrast,  uses  con¬ 
cepts  from  genetic  algorithm  search  routines  to  more  efficiently 
search  this  parameter  space.  The  sampling  is  done  in  such  a  way 
as  to  avoid  becoming  stuck  in  locally  optimal  solutions.  The  effi¬ 
cacy  of  the  approach  has  been  demonstrated  in  identifying  crack 
location,  length,  and  orientation  using  only  simulated  impulse-re¬ 
sponse  data  from  a  cracked-plate  model.  Because  the  Pop-based 
MCMC  method  requires  repeated  iterations  of  the  forward  model, 
great  care  was  taken  in  developing  an  efficient  model.  This  was 
accomplished  by  using  tailored  “serendipity”  elements  to  describe 
the  stress  field  near  the  crack  tip.  This  allows  for  good  model  con¬ 
vergence  with  many  fewer  elements  than  are  used  in  standard  fi¬ 
nite  element  codes.  The  combination  of  efficient  modeling  and 
effective  parameter  identification  routines  can  provide  a  wealth 
of  information  about  the  state  of  a  structure.  Both  the  parameter 
estimates  and  the  credible  intervals  associated  with  those  esti¬ 
mates  are  obtained.  This  allows  for  confidence-based  decisions 
regarding  the  maintenance  of  a  structure  and  also  provides  infor¬ 
mation  needed  in  prognostics  models  for  damage  evolution  in 
structures. 
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